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1 . GENERAL 

Recently, non-axisymmetric convection in (vertical) directional solidification 
experiments has been observed. It has been suggested that the flow character is a 
consequence of the lack of azimuthal symmetry in the temperature field. Motivated by 
these observations we have examined the consequences of deviations from axisymmetric 
wall temperature conditions in a vertical differentially heated cylinder. We show that the 
degree of asymmetry exhibited by the flow depends on the ratio between the amplitude of 
the maximum azimuthal and vertical temperature difference and on the ratio between the 
maximum calculated fluid velocity magnitude and “thermal diffusion’ speed (thermal Peclet 
number). 

The object of this work was to develop and implement a numerical method capable of 
solving the non-linear partial differential equations governing heat, mass and momentum 
transfer in a 3D cylindrical geometry in order to examine the character of convection in an 
asymmetrically heated cylindrical ampoule. The details of the numerical method, including 
verification tests involving comparison with results obtained from other methods, is given 
in Appendix 1. The results 1 of our study of 3D convection in an asymmetrically heated 
cylinder is described in the following sections. 

2. INTRODUCTION AN BACKGROUND 

The character of convection during solidification has been examined using 
numerical models of buoyancy-driven convection in cylindrical and rectangular geometries 
[1-12]. Early work [1-3] includes a variety of imposed temperature boundary conditions. 
These range from purely vertical temperature gradients which results in convection after a 
critical value of the Rayleigh number is exceeded [4], to idealized conditions associated 
with Bridgman-Stockbarger furnaces [3] which are imposed directly on the melt and crystal 
without consideration of the heat transfer between the ampoule, furnace and sample. For 
these boundary conditions, flow always occurs owing to the presence of radial temperature 
gradients. Later models [5,6] have accounted for the presence of the ampoule, and the 
details of furnace design. In an actual growth situation the thermal profile of the inner 
surface of the furnace is not realized at the ampoule wall; it is modified by heat transfer 
between the crystal, melt, ampoule and the furnace itself. The tendency is to reduce axial 
temperature gradients, while radial temperature gradients may increase or decrease 


lr These results are contained in paper submitted to the Physics of Fluids A. 
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depending on the specific nature of the heat transfer between the charge and ampoule [5]. 
The influence of melt convection on dopant distribution has been examined for dilute and 
non-dilute melts [3,5,9]. For a given furnace- ampoule combination the amount of 
compositional non-uniformity (or radial segregation) was shown to be a non-linear function 
of the Rayleigh number “Ra”. The growth rate and physical properties of the melt also 
influence the degree of compositional uniformity in the grown crystal. These studies were 
restricted to cases for which the thermal boundary conditions possessed azimuthal 
symmetry. 

Recently, non-axisymmetric convection has been observed in directional 
solidification experiments conducted with low melting point, low thermal conductivity 
materials [13,14], The observed flows were explained by lack of azimuthal symmetry in 
the temperature field caused by the fact that the particular Bridgman-Stockbarger 
solidification apparatus precluded any control of the azimuthal temperature field. Indeed, 
such a situation may be common in many instances where the Bridgman-Stockbarger 
technique is employed owing to the difficulties involved with constructing an axisymmetric 
heating arrangement and in aligning the ampoule axis with the axis of the furnace. During 
growth of high temperature materials, any misalignment of the ampoule will result in an 
azimuthal variation in the radiation view factors for the ampoule. This will lead to 
asymmetry in the temperature of the ampoule wall. 

The object of this work is to examine the character of convection in an 
asymmetrically heated cylindrical ampoule. In section 3 we formulate the model problem. 
In section 4 we outline the pseudo-spectral collocation method used to solve the 3D 
problem. The results are presented in section 5 and discussed in section 6. 

3. FORMULATION OF THE PROBLEM 

A practically desirable Bridgman-Stockbarger set-up such as that described by 
Dahkoul et al. [15] consists of three distinct thermal zones. The simplest arrangement has a 
hot zone (in which, ideally, the temperature at the wall and in the melt should be almost 
isothermal) and a cold zone. These are separated by a gradient zone in which a thermal 
barrier controls the heat transfer between the ampoule and the hot and cold zones. The basic 
Bridgman-Stockbarger model used here is shown in Fig. 1 . The fluid is differentially heated 
in a vertical cylinder having a length 2H and a diameter D. The lower and upper endwalls 
of the cylinder are maintained at constant temperatures Tm and Th respectively (where Tm 
< Th)- In this system the presence of the thermal barrier is approximated by taking the 
cylinder walls to be adiabatic in the region corresponding to the gradient zone. In the hot 
zone we impose asymmetry in the temperature distribution using a function which causes 
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azimuthal and vertical deviations from isothermal conditions on the cylinder wall. The 
maximum lateral deviation in temperatue is given by ATq (see Fig. 1). All boundaries of the 

cylinder are rigid with no-slip conditions for the velocities. 

The equations governing energy, momentum and mass transport are written in 
primitive variables for cylindrical coordinates and using the Boussinesq approximation . 
They are rendered dimensionless using L*, V*, L*/V* and Th - Tm as characteristic 
scales for length, velocity, time and temperature, respectively. The 
dimensionless equations are: 


— = — V • grad T 

9t Re Pr 


— = .M. . v) . y - grad p + — T k. 

Re Pr Re 2 

div V = 0, 


( 1 ) 

( 2 ) 

( 3 ) 


where t is the time and V = ( v r , v 0 , v z ), p and T =^^ respectively represent the 
dimensionless velocity, presure and temperature. The unit vector in the z-direction is 


denoted by k, Pr = is the Prandtl number, with v the kinematic viscosity and K the 
K 

thermal diffusivity, Ra = ga(L * ) ^J H ' T ^ is the Rayleigh number, with g the gravitational 

acceleration, and a the thermal expansion coefficient, and Re - L*V* i s the Reynolds 
number. 

Several choices are possible for these characteristic scales. In our computations, 

after trying several kinds of scaling, the choice V* = (i.e. Re = ^ = Grashof 

number) proved to be the best for obtaining good convergence behavior. This scaling has 
been suggested by Ostrach [16]. For the choice of L*, the length H proved to be more 
convenient for computational purposes, primarily because of the chosen numerical method 
(see section 3). By exploring the different scalings, we have noted that the convergence 
behavior of the numerical method described in section 3, could be changed significantly 
(for the worse) when the parameters multiplying diffusive terms (Laplacians) are too small 
or too large. These considerations led us to recast equations (1) and (2) in the form 
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— = -LAT-|a-V-grad T , (4) 

dt Pr Pr ~ 

— - = AY - ^4- (grad V)* V - grad p + T k. (5) 

at Pr 

Finally equations (3-5) are solved as functions of the independent variables (r,0,z) in the 
domain Q = ] 0, -T [ x [ 0, 2 jc [ x ] -1, 1 [. Here the aspect ratio A is defined as 2H/D, 

A 

where D is the diameter. 

The following conditions were applied at the initial iterate (t < 0) 

v r = ve = v z = 0, 

T = g(r,6,z) f(z), 

wi,h f<z)= i^^[ exp (^)' exp(a) 

and g(r,0,z) =1 + A 0 ^ expj^l) (cos(7t-6)-l) sin^tJ^A). 

In order to avoid discontinuities in derivatives of the temperature, the abrupt transition 
between the hot zone and the adiabatic zone the function f(z;a) is used. The parameter a is 
chosen such that were it not for the azimuthal temperature perturbation, the dimensionless 
wall temperature due to f(z;o) would be almost equal to 1 for z > Z a where Za is the z- 
coordinate of the location of the upper limit of the adiabatic zone. The function g(r,9,z) 
gives the initial condition for the azimuthal variation in temperature. The maximum 
azimuthal temperature deviation is denoted by ATq. It is convenient to express this 

rr^ rr^ • _ * AT0XlOO% 

deviation as a percentage of the axial temperature difference Tr - Tm; i.e. Ae = Tr . t m ' 

Table 1 gives Ae (A e =2 Aoj^jexpj^Aj ) for different values of A 0 and different aspect 
ratios. 

The following dimensionless boundary conditions are applied at t > 0 


v r = ve = v z = 0 at z = ± 1 and r = Jk 

(7a) 

T = 0 at z = -1, T = 1 at z = 1, 

(7b) 

T = g(r,9,z)f(z) at r = A- and z > Z a (z ±1) 

r\ 

(7c) 


(6a) 

(6b) 
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2i=0at r = l andz<Z a , C7d) 

dr A 

The azimuthal temperature variation on the ampoule wall is given by g(j^,0,z) and is 
enforced only for the hot zone z > Za- 

Table 1. Values of the percentage of azimuthal temperature difference Aq 


Aspect Ratio 

> 

o 

Ae 

1 

0.1 

10% 

1 

0.02 

2% 

1 

0.2 

20% 

2 

0.1 

13.2% 

2 

0.2 

26.4% 


4. NUMERICAL METHOD 

The equations (3-5) have been solved using a modified version of the Fourier- 
Chebyshev pseudospectral method introduced by Pulicani and Ouazzani [17]. For this 
method the singularity which arises at r=0 when using cylindrical coordinates (without 
axisymmetry) [18] is avoided by using a change of dependent variables. The singularity 
arises because r=0 is an artificial boundary of the computational domain; i.e. it occurs only 
by construction of the coordinate system. The change of variables used to cope with the 
singularity is 


T = 


T 

r 



and 



( 8 ) 


Application of (3-8) yields a system of equations with T, v T , v^, p and ve, as unknowns. 
Obviously, T = v^ = v^ = v^ = p = 0atr = 0. These equations are then discretized with 
respect to time, t = n 8t, by means of a second order semi-implicit scheme. The latter 
employs a combination of the Adams-Bashforth and Crank-Nicolson schemes, namely 
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with 



(9) 


r .L|lA 


?n+l 


- v r 


~n,n-l,n+l n 

= G - 5t r [grad j 


n+1 


( 10 ) 



( 11 ) 



= 5t 



_ 3 R a y 
2 Pr 




and 


~n,n-l,n+l [ 

G = St r A 


-^-(gradyl-V 


||a[(gra<l|)-V ] n l + [Tk]""', 


where for convenience we have defined V = ( v r , ve, v z ). 

With this method [17], problems encountered with satisfying (11) are surmounted 
by using artificial compressibility [19,20], A false timestep is employed and the method is 
only applicable for steady solutions to the system (3-7). Because of the stiffness of the 
physical problem at low Prandtl numbers, obtaining rapid convergence with this method at 
high Rayleigh numbers is difficult. In order to overcome this problem, we have modified 
the method by introducing two iterative processes. An outer iteration which is related to 
each timestep and an inner iteration which ensures that (11) is satisfied to some e«0 for 
each outer iterative step. It is clear that for lower values of e, this pseudo-unsteady method 
takes on the character of an unsteady calculation. For this reason we have used an Adams- 
Bashforth scheme to discretize the convective terms instead of simply taking them explicitly 

~n,n-l ~n,n-l,n-Fl 

[17] at the instant t = n St in F and G 

A generalized ADI procedure [19] is then applied to reduce the problem to the 
successive solution of one-dimensional problems. For clarity, we present the method only 
as applied to the momentum transport equation. It is readily extended to cope with the 
energy transport equation. At each time step the following problem is solved 



A 0 )v*=G n ’ n ' 1,n+1 


5tr 


[grad j 


n+l,|l 


(12a) 
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(l-|A,)v. V*. 

(12b) 

(r-tflA,)(r>V*) -V”. 

(12c) 

y n+l.ji+1 _ y n + y***^ 

(12d) 

~n+ 1 ,n.+t _ ~n+t,n _ ^ r ( div Y ) n+l4l+1 5 

(12e) 


P 

grad . Equation 


where X is a strictly positive constant and p (p = 0,1,. ..Np) is the superscript 
connected to the inner iterative process. The operator A occurring in (10) has been 
decomposed so that A = A r + Aq+ A z . Here the A®, are related to the independent 

r p"ln+l,0 

variables r, 0 and z, respectively. Note that for p = 0, [grad — 

(12e) has no physical meaning until 5p n+1 = p n+1,R+1 - p n+1 -n = 0 (or practically 
speaking, is sufficiently small). If 5p n+1 is identically zero at each time step the method is 
no longer pseudo-unsteady. The present way to deal with the pressure can be time 
consuming when the number of internal iterations Np is too high. However, since we seek 
steady solutions, the iterative process has been introduced only to help the convergence 
when the solution is very stiff. It is stopped as soon as the divergence of the velocity 
reaches a certain value £. The maximum value of Np is chosen such that the convergence is 
obtained after a reasonable number of outer iterations, but not so high that the divergence of 
the velocity is £ at the beginning of the outer iteration. The optimum values of Np and £ will 
be defined in the next section. A standard Fourier-Galerkin approximation [18] is 
employed for the solution of (12a). 

The most convenient (see section 2) choice of characteristic length is L* = H. This 
leads directly to the domain -1 <, z < 1 which is required for the Chebyshev-collocation 
method [20] used to solve equations (12b) and (12c). The boundary conditions are 
introduced by replacing the right hand sides of (12b) and (12c) by the appropriate terms 


[17]. 

The energy equation is solved in the same manner as the momentum transport 
equations but, obviously, needs no inner iterations to satisfy the divergence equation. The 
solution algorithm takes the following steps : 

(i) With T n , T n_1 , V n and V n_1 known, we deduce r and then T n+1 . 

(ii) Using T n+1 , V n and V n '\ we calculate G"'" ’" + and then, finally, V n+1 and 
p n+1 with (12). 
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In all of the results we present here the above method requires a distribution of 
points such that 


e k =^-, k = o,..,N e -i, 


Zj = cos|jH_), j = 0,..,N Z -1, 

ri= ^Hsw) +1 . 


, i = 0,..,N r -l, 


(13a) 

(13b) 

(13c) 


with N r , Ng and N z the number of collocation points between r = 0 and 1/A, 0 = 0 and 
2 tc, and z = -1 and 1, respectively. 

In our discussion of the results we shall also refer to the residual of ft, denoted R^, 
which is calculated on the collocation points of ft = (T, v r , v a , v z }. Furthermore, 


R d = Max iik j 


1 'fl(r i ,9i c ,z i ) n+1 - t3(ri,9 k , Zj 
St t3(ri,0 k ,Zj)r +l 
with i = 0,..,N r -l * k = 0,..,N k -l > j = 0,..,N Z -1 


(14) 


We denote the maximum values of v r , vq, v z , v x and v y by v r max , vg max , v z max , v x max 
and v y max , respectively. The variables v x and v y are the velocities in Cartesian coordinates 
such that v x = v r cos 0 - vg sin 0, and v y = v r sin 0 + vg cos 0. As these maxima are 
calculated at the collocation points particular to each method, small differences are expected 
in the results. 

For all the results presented in the following section the starting condition (6b) has 
only been used for Ra = 10. For the others, the starting condition was the solution 
calculated with a lower Ra. The adiabatic zone covers one fourth of the cylinder's height 
when the aspect ratio A = 1 (i.e. Z a = -0.5) and one eighth when A = 2 (i.e. Z a = -0.75). 
The calculations have been performed on a Cray XMP computer by using a spatial 
resolution such that N r x Ng x N z = 13 x 20 x 31 when A = 1 and N r x Ng x N z = 13 x 
20 x 61 when A = 2. For all calculations the parameter X (see (12e)) is equal to 1.3, 20< 
Njx< 30, and the time-step 8t chosen between 1. 10" 4 and 1. 10"^ according to the stiffness 
of the solution. Note that we have stopped our calculations when the divergence div V « 
10' 8 and ResV = 10" 4 . No significant change in the values of the velocity and temperature 
fields if divV and ResV are decreased further. For the problem described here a 
comparison has been made between our method and the finite element code FIDAP 
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[21,22]. For FIDAP an irregularly spaced Cartesian grid with nodes N x xN y xN z = (9 x 
llx 31) was employed (see Fig. 2). 

For the purposes of presenting the results we define two dimensionless measures of the 
velocity, a Reynolds number Re* ma x = Re*max = v* max D/v and a thermal Peclet number 
Pe* m ax = v m axD/K 7 w here a corresponds to the coordinate directions (i.e. r, z, q, x, or 
y). The thermal Peclet number represents a ratio of the magnitudes of the maximum 

dimesional convective velocity and a thermal diffusion velocity. 

5. RESULTS 

Table 2 shows the results of a comparison between our method and the finite 
element code FIDAP [21,22] . The Values of Re x max , Re y max an d ^ e zmax> were obtained 
with Ra = 250, 2500 and 1500. Clearly, there is good agreement between these two 
methods for Ra between 250 and 2500, while for Ra= 15000, the poor agreement is only 
due to the smaller number of points in the z-direction employed for FIDAP. 

Table 2. Comparison between results obtained from the Spectral method and from FIDAP 



SPECTRA1 

L 

FIDAP 


Ra 

Re x max 

R^vmax 

Re z max 

Re x max 

Re v max 

Re z max 

250 

10 

7 

20 

12 

8 

21 

2500 

60 

61 

103 

067 

069 

104 

15000 

263 

192 

286 

295 

216 

386 


Table 3 summarizes the details of our computations for A — 1,2, 2500 < Ra ^ 64000 and 
2% < A 9 < 27%. A comparison of the velocity fields in the 0 = 0° and 180° sections for 
A 0 = 2% and Ra=2500, and Ra =24,000 is given in Fig. 3. The flow is barely perceptible 
at the lower Ra. Comparison of Fig. 3 with Fig. 4., which shows the velocity fields for Ra 
= 2500, 15000 and 24000, reveals the effect of increasing the temperature asymmetry. The 
effect of increasing Ra at fixed Ag is also seen in Fig. 4. Note that the locations of the roll 
centers change as Ra is increased. Figure 5. depicts the velocity and temperature fields for 
three vertcial sections at Ra =15,000 for A 9 = 20%. Horizontal sections of the velocity field 
for the Ra= 15,000 case are shown in Fig. 6 . 

The effect of increasing the aspect ratio is seen upon comparison of Fig. 4 with 
Fig. 7 which has been calculated for Ra = 2500 and 15000 with A 0 = 0.2 (Ag = 26.4%). 
The basic asymmetry' of the flow is not affected significantly by the increase in aspect ratio, 
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although the centers of the toroidal rolls are shifted. Horizontal sections of the flow shown 
in Fig. 6 are depicted in Fig.7. 

Table 3. Summary of the results for Pr = 10" 2 


Ra 

A 

> 

o 

5t 

max 

R^Gmax 

max 

2500 

1 

0 

10- 4 

70 

0 

122 


1 

0.02 

10- 4 

69 

7 

120 


1 

0.2 

10- 4 

61 

56 

103 


2 

0.1 

10-5 

69 

44 

114 


2 

0.2 

10-5 

69 

59 

125 

6400 

1 

0.2 

T 

o 

112 

123 

168 


2 

0.1 

10- 4 

128 

99 

191 

15000 

1 

0 

10- 4 

205 

0 

314 


1 

0.02 

10- 4 

222 

42 

307 


1 

0.2 

10- 4 

220 

254 

286 


2 

0.1 

lO- 5 

209 

175 

294 


2 

0.2 

10*5 

188 

228 

425 

24000 

1 

0 

10- 4 

242 

0 

384 


1 

0.02 

10- 4 

267 

62 

376 


1 

0.2 

7x10-5 

318 

363 

373 

64000 

1 

0 

10.- 4 

338 

0 

583 


Comparison of Figs. 8 and 9 for Ra = 15000 and Ae = 13.2% with Figs. 10 and 1 1 
for Ra=2500, Ae =13.2%, reveals that in the higher Rayleigh number case two additional 
cells have formed in the upper half of the cylinder. These are barely detectable in the 
Ra=2500 case. At higher values of Ae these cells are able to develop well for lower values 
of Ra. In Fig. 7 the twofold increase in Ae leads to well developed upper cells even at Ra 
= 2500. 
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Figs. 12 and 13. depict the results obtained for Pr = 1, Ra = Gr = 2.5x10 s . 
Comparison of these results with Figs. 8 and 9 (Pr = 10'2, Ra = 2500, Gr = 2.5x10^) 
reveals the effect of increasing Pr while holding the Grashof number fixed. There is a 
decrease in the degree of asymmetry in the flow for the lower Prandtl number case. The 
cylinder is dominated by a large asymmetric roll which extends along most of the cylinder 
with smaller secondary roll along the bottom of the cylinder. The isotherms which extend 
into the adiabatic zone have been modified by the flow and are considerably flatter than 
their low Pr counterparts. 

6. DISCUSSION AND SUMMARY 

The results of an experimental investigation into the nature of asymmetric flow 
during Bridgman-Stockbarger directional solidification of Salol by Neugebauer and Wilcox 
[14] led to the conclusion that the degree of flow asymmetry decreases with increasing 
convective flow velocities. Furthermore they conjectured that for the same experimental 
conditions, low Prandtl number fluids to exhibit correspondingly less asymmetry owing to 
the more rapid flow velocities which would necessarily ensue. It can be discerned from the 
results presented in section 4 that our calculations predict that for a fixed value of A 0 and Pr 
an increase in Re* max (in practical terms caused by an increase in AT=Tpj-Tm ) tends to 
amplify the asymmetry. Examination of the experimental results of [14] reveals that when 
AT was increased, ATq, the azimuthal variation in temperature, remained the same. In other 
words, the relative temperature asymmetry Aq was decreased. Thus, the observed decrease 
in flow asymmetry can be explained merely by the fact that the azimuthal temperature 
variation was less significant for the higher Ra cases. This trend is confirmed by our 
calculations. 

Our calculations do confirm the prediction made in [14] that, for otherwise 
equivalent conditions, the higher Pr case will exhibit more asymmetry in the flow than the 
low Pr case. That the flow asymmetry is reduced due to an increase in flow velocity 
(Reynolds number)contradicts our results presented in Table 3 and Figs. 4-9. The 
reduction in asymmetry can be explained, however, upon examination of relationship 
between the degree of flow asymmetry and the thermal Peclet number. The latter scale 
represents a thermal velocity magnitude. Table 4 clearly shows that at a fixed Aq there is an 
increase in the degree of flow asymmetry as the thermal Peclet number increases. 
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Table 4. Comparison of Peclet numbers for the Pr = 10' 2 and Pr =1 cases. 


Pr 

P e r max 

P e 0 max 

max 

lO- 2 

0.69 

0.44 

1.14 

1 

11.02 

8.92 

22.83 


In summary, we have examined the consequences of azimuthal asymmetry in 
ampoule wall temperature in a differentially heated cylindrical ampoule containing an 
incompressible Newtonian fluid. The study was motivated by recent observations of non- 
axisymmetric convection in directional solidification experiments conducted with low 
melting point, low thermal conductivity materials [13,14]. Our results indicate that for a 
fixed value of the relative asymmetry in temperature Ae, the degree of asymmetry in the 
flow is accentuated as the thermal Peclet number is increased. For fixed Pr and Ra, the 
flow asymmetry increases if Aq is increased. 
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8 . FIGURE CAPTIONS 

Fig. 1 Model Bridgman-Stockbarger configuration. 

Fig. 2 Mesh used for the FIDAP calculations. 

Fig. 3 Comparison of the velocity fields in the 0 = 0°and 180° sections and the (r,z) 
plane at z = 0 with A 0 =2% and Ra-2500, 24000. Aspect ratio A=1 and Pr=10- 2 . 
Fig. 4 Comparison of the velocity fields in the 0 = 0°and 180° section with A 0 =2O% 

and Ra=2500, 15000, and 24000 . Aspect ratio A=1 and Pr=10' 2 . 

Fig. 5 Velocity and temperature fields with A 0 =2O% and Ra=15000 for 

0=0° and 180, 54° and 234°, 90° and 270°. Aspect ratio A=1 and Pr=10- 2 . 

Fig. 6 Velocity field in the (r, 0) plane with A 0 =2O% and Ra=15000 at z=0.75, 0, 

and -0.75. Aspect ratio A=1 and Pr=10' 2 . 

Pig # 7 Velocity fields with A 0 = 26.4% Ra = 2500 and 15000 at 0 = 0 and 180, 54 and 

234°, 90° and 270°. Aspect ratio A=2 and Pr=10‘ 2 . 

Fig. 8 Velocity and temperature fields with A 0 =13.2% and Ra= 15000 at 

0= 0° and 180, 54° and 234°, 90° and 270°. Aspect ratio A=2 and Pr=10‘ 2 . 

Fig. 9 Velocity field in the (r, 0) plane with Ra=15000 at z=0.5, 0, -0.5, and -0.87. 
Aspect ratio A=2 and Pr=10‘ 2 . 

Fig. 10 Velocity field with A e =13.2% and Ra=2500 at 0=0°, 54°, and 90°. Aspect ratio 
A=2 and Pr=10' 2 . 

Fig. 1 1 Velocity field in the (r, 0) plane with Ra=2500 at z=0.5, 0, -0.5, and -0.87. 
Aspect ratio A=2 and Pr=10' 2 . 

Fig. 12 Velocity and temperature fields with A 0 =13.2% and Ra=250000 at 

0=0° and 180, 54° and 234°, 90° and 270°. Aspect ratio A=2 and Pr=l. 

Fig. 13 Velocity field in the (r, 0) plane with Ra=250000 at z=0.5, 0, -0.5, and -0.87. 
Aspect ratio A=2 and Pr=l. 
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9. APPENDIX 1: 

A Fourier-Chebyshev Pseudo-Spectral method for Solving Steady 3D 
Navier-Stokes and Heat Equations in Cylindrical Cavities, by J. P. Pulicani 
and J. Ouazzani (submitted to Computers and Fluids, August 1990). 
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A Fourier-Chebyshev Pseudospectral Method for Solving Steady 3-D 
Navier-Stokes and Heat Equations in Cylindrical Cavities 

By 

J.P. Pulicani and J. Ouazzani 

Center for Microgravity and Materials Research, University of Alabama in Huntsville, 

Huntsville, Alabama 35899, USA 

Abstract 

A Fourier-Chebyshev pseudospectral method for solving the steady 3-D Navier-Stokes 
and energy equations in cylindrical cavities is presented and discussed. The general method is 
pseudo-unsteady and uses a semi-implicit finite difference scheme for the time integration. The 
generalized ADI procedure is then applied to reduce the problem to successive solutions of 
one-dimensional problems. The spatial discretization uses a Fourier-Galerkin approximation in the 
periodic direction and a Chebyshev-collocation approximation in the other directions. Difficulties 
related to the pressure are surmounted by using the artificial compressibility method. A suitable 
variable change has been chosen to avoid the problem of singularity at the axis generated by 
cylindrical coordinates. The method is first tested on an advection-diffusion equation and then on 
the Navier-Stokes and heat equations. Finally, the method is illustrated by the problem of 
convection in a differentially heated fluid. 
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1. Introduction 

For a long time numerical fluid mechanics has been focused mostly on two-dimensional 
models. Limitations related to the computer performance (i.e. excessive CPU time and large 
memory storage) have been among the main reasons for this trend, which has often prevented the 
solution of three-dimensional problems. However, two-dimensionality is difficult to achieve in real 
experiments. Nevertheless, increased research on numerical analysis and the improvement of 
computer performance permit more efficient solution of the complete 3-D Navier-Stokes and 
energy equations for several cases. Of the many numerical methods in use, those based on spectral 
approximations have met with increasing success over the last ten years. The principal attributes of 
spectral methods (see [1-4]) are their high accuracy and dual representation of the dependent 
variables in both physical and spectral space. Another advantage of these methods is their ability to 
accurately represent the solution with fewer collocation points than others (i.e. finite difference, 
finite elements, etc). This last advantage makes 3-D problems more tractable when using spectral 
methods. However, even though these methods have been extended to a large variety of 2-D 
physical problems (such as turbulence, channel flows, natural convection, crystal growth 
processes, etc) they have rarely been used for 3-D problems involving a cylindrical geometry. 
Most models assume axisymmetric conditions which result in 2-D problems. 

In this paper we present and discuss a Fourier-Chebyshev pseudospectral method to 
solve the incompressible 3-D Navier-Stokes and energy equations. These equations are solved in a 
cylindrical geometry without the assumption of axisymmetry. The general method is 
pseudo-unsteady and uses a semi-implicit finite difference scheme for the time integration. The 
generalized ADI ( 'mating Direction Implicit) procedure [2] is then applied to reduce the problem 
to the successiv .uon of one-dimensional problems. The spatial discretization uses a 
Fourier-Galerkin approximation in the periodic direction and a Chebyshev-collocation 
approximation in the other directions. Similar techniques were introduced for the solution of 2-D 
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steady binary gas mixture flows [5], and for the solution of a diffusion equation in a domain 
between two concentric spheres [3]. Difficulties related to the incompressibility condition and the 
pressure are surmounted by using the artificial compressibility method as used in [2,5]. Another 
method to maintain the incompressibility constraint is the use of the influence matrix technique [6]. 
This latter technique was applied for incompressible flows in cartesian and cylindrical geometries 
by [7] and in a rotating annulus by [8]. When using cylindrical coordinates without axisymmetry, 
the complexity is further increased by the presence of a singularity at the axis [7]. To eliminate this 
singularity, various methods such as the use of Cartesian coordinates [9,10], variable change [11], 
iterative technique [12] and special mesh [13], have been applied. In this paper we propose and 
assess two kinds of treatment. The first one is based on a variable change and the other on a device 
to avoid the points at the axis. 

The method is first described and tested on an advection-diffusion equation and then on 
the Navier-Stokes and heat equations. Finally, the method is illustrated by a convection problem of 
a fluid in a horizontal differentially heated cylinder. 


2. Numerical method for an advection-diffusion type equation 


2.1 General method 


For the description of the method, we consider the following problem : 


— = AT - V.VT + F t 
3t 


, t>0, 


(2.1a) 


4 


where, 


AT = 


d 2 ! 

dr 2 


+ 


1 * + 

r 9r 


t d 2 T | 3 2 T 
r2 30 2 9z 2 ’ 


(2.1b) 


V.VT = v 



r 


3T 

— + v z 

ae 


3T 

8z 


V = ( v r , ve, v z ) , 


(2.1c) 


for (r,9,z) e Q = ]0, R[ x [0, 2ji[x ]— 1, 1[ (see Fig. 1) and the boundary and initial conditions : 

ap Tip + bp — I = cp , on F = 3f2, (2.2a) 

3n T 

T = T 0 , at t = 0, (2-2b) 

where ap, bp, cp, T 0) v r , vq, v z and Fj are given functions and d /9n is the normal derivative. In 
real physical problems, V is the velocity. 

Equation (2.1a) is discretized with respect to the time, t = n 8t, using a semi-implicit 
scheme for the diffusive terms and explicit for the convective term : 

(l-a5tA)(T n+1 -T n ) = 5t F 1 , (o>0) , (2.3a) 

= ( AT - V.VT + F t ) n . (2.3b) 

Note that a = 1/2 corresponds to the Crank-Nicolson scheme. 

The discrete equation (2.3) is solved by means of the generalized ADI procedure (as used in [2,5]) 
which reduces it to a set of one-dimensional problems. At each time step the following problem 
will be solved : 
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( l - o 5t Ae ) 'P 

= 5t F" , 

(2.4a) 

s|c 

( l-aStAz)'? 

* 

= V P , 

(2.4b) 

*** 

( 1 - a 5t A r ) *F 

** 

= V P , 

(2.4c) 

q-n+l 

*** 

= 1^ + ^ , 

(2.4d) 

Ae=> V Az = 

8 , 8 i 8 

, 2 ’ r %2 + r 9r ' 

(2.4e) 


A standard Fourier-Galerkin approximation is employed for the solution of (2.4a). Thus, 
'P is expanded in a truncated Fourier series as follows 

V*(r. 9, z) = ^ 'FkCriz) e* 6 , F n (r, 0, z) = ^ Fg(r,z) e lke . (2.5) 

k=-K/2 k=-K/2 

Equation (2.4a) can then be written as 

'Pk(r,z) + a 8t^-'Fk(r,z) = 5tFk(r,z) , (for each harmonic k) , (2.6) 

r 2 

with the Fourier coefficients 'Fk as unknowns. From the solution of (2.6), it is then 
straightforward to deduce with (2.5). The values of V P* have been calculated at the collocation 
points i.e. 

9k = 27tk/K , k = 0,...,K-1. (2.7) 

The equations (2.4b) and (2.4c) are solved by means of a Chebyshev-collocation method. 

For 'P , we apply a Chebyshev polynomial expansion in the z-direction : 
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M 


Y (r,9,z) = X (r,0) T m (z). 

m=0 


( 2 . 8 ) 


where refers t0 the m th Chebyshev coefficient and T m to the m th Chebyshev polynomial. 
The orthogonal collocation consists in expressing the derivatives at one collocation point in terms 
of values of the function at all points. For instance the second derivative is expressed by 



at the collocation points : zj = cos(7tj / M) , j = 0,...,M. (2. 10) 

The coefficients d(zj,z^)( 2 ) are obtained by using the orthogonality property of the Chebyshev 
polynomials and the usual trigonometric formulas [3,5]. In matrix-vector notation, the equation 
(2.9) can be written as 

2 ** 

— ^r-(r,0,Zj) = , © = [ d(z j ,z JX )( 2 >] , (2.11a) 

dz 2 

0 = [ x F**(r,0,zo) (r.G.Zp.) (r,0,ZM) ] T - (2.11b) 


Then by using (2.11) in the left-hand side of (2.4b), we obtain the algebraic system 
$ # ♦ 

<|> = <[) , 5Vf = I - cr 8t T> , (2. 12a) 

<t>* = [ V(r,0,zo) , ... , ^(r.e.zn) , ... , ^*(r,0,z M ) f- ( 2 - 12b ) 

For the solution of (2.4c), the technique is identical, but a Chebyshev polynomial expansion is 

*** 

applied in the r-direction for 'F , that is 
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N 


WWW *** 

V (SM = X (0,Z) T n (^) , 

n=0 


(2.13) 


where £, refers to the transformed plane ^ = 2r / R - 1 (-1<^<1 and 0 < r < R ). 
The first and second derivatives of *P are then written as 


iPw N , » **» 

<^i, 0 ,z) = £ a(^) (p) *P (%n.0.z) , P= 1,2 
n=o 


(2.14) 


5i = cos(7ri/N) , ri = R (^i - 1) / 2, i = 0.....N , 


(2.15) 


( <•) *** ***^ 

i a 1 ? „ ft , 

■ - + ^~ - — l(ri,0,z) = 


LaV 


ar 


2 ri 0 r 


*** *** 1 


+ — (^i,0,z) = * 0*** , (2.16a) 

-\e 2 Ti -.e 


af 11 ^ 


Jl = ». 2 [a&^) <2> ] + ^[a«i,^) (1> ] . * = (r) 


(2.16b) 


*** *** 


0 =[ x p (^o.0.z) 'p (^.0.z) . - . 'p (^N>e,z) r. 


(2.16c) 


We obtain an algebraic system by using (2.16) in the left-hand side of (2.4c) in the form : 

Jjtsjc 

^0 =0 , 5Y> I - o St ^ , (2.17a) 

♦” = [ >p**(4o.e,2) , ... , <p**«,,e.z) >p'*ft N ,8,z) ] T . (2.17b) 

The first and last equations in (2.12) and (2.17) are replaced by the boundary conditions. If the 
boundary conditions are not time-dependent and if the initial condition satisfies these boundary 
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** *** 

conditions, then the boundary conditions for <|> and <t> are homogeneous. Still that in all cases, 
the inversion of the resulting matrix M and is made only once before the time integration 

begins. Thus, only matrix vector products have to be performed at each time-step, that is : 

** * *** ** 

<t> = <(> , then <(> = 9t l 4> and finally T n+1 with (2.4d). 

The right hand side of (2.4a) is evaluated by the standard pseudospectral technique [1]. The 
derivatives are performed in the spectral space and the products in the physical one. These two 
spaces are connected through a EFT algorithm [14]. Note that the FFT algorithm is more efficient 
than the direct matrix-vector multiplication method only when the number of collocation points is 
greater than 30. 

When using cylindrical coordinates there is a difficulty associated with the presence of 
the singularity at the axis (for 1/r — >°° when r — >0). This difficulty arises because the axis is a 
boundary of the computational domain only by construction of the coordinate system. Obviously, 
this problem can be alleviated when the value of T(r*=O,0,z) or its derivative is known at the axis 
(for instance if the solution is supposed to be axisymmetric or if a natural boundary condition 
exists to be imposed at this axis). For most physical problems, the axis cannot be considered as a 
boundary in the problem (2.1); therefore a special treatment to eliminate this singularity is required. 
In this paper we propose two kinds of treatment. One is based on a variable change and the other 
on a distribution of points exluding the axis. 

2.2 Treatment of the singularity by using a variable change 

The problem of singularity can be avoided by applying to equation (2.1) the following 
change of variable : 


T(r,0,z) = rP T(r,0,z) , (3 < 0. 


(2.18) 
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The idea for such a variable change came from a private conversation with Dr R Peyret (University 
of Nice, France). A similar variable change for a functional equation has also been used in [15, 
p. 1104]. The main advantage of the variable change (2.18) is that T(r,0,z) = 0 at r = 0 for P < 0. 
For the problem (2.1)-(2.2) it is normally not necessary to use the same variable change for the 
velocities since they are given functions. But in order to apply the method to the Navier-Stokes and 
energy equations, we enforce the following identities : 

v r (r,0,z) = rP v r (r,0,z), v e (r,0,z) = rP ^(r,0,z) and v z (r,0,z) = rP v z (r,0,z). (2. 19) 


Applying these variable changes to (2.1a) it yields 


— = F = r - P ( AT - V.VT + F T ) , (2.20a) 

0t 


r 


-P AT = 


3 2 T (2p+i) af 
ar 2+ r * + 



i a 2 T | a 2 T 

r2 ao 2 az 2 


(2.20b) 


r-P V.VT = rP v T — - + P rP' 1 v x T + rP- 1 v 0 — + rP v^ — 
9r 30 dz 


(2.20c) 


Thus, equation (2.20) remains to be solved for T. Furthermore, due to the change of variable 
(2.18), the boundary and initial conditions are now : 


Tl 4. h 

ari T 'r, + b ri U p — -T + — 
\ On 9n 


= r-P 


c n 


on Ti 


(2.21a) 


Tlr 2 = 0 , on f 2 


(2.21b) 


10 


T = r-PT 0 , at t = 0 


(2.21c) 


where I ~2 refers to the boundary r = 0 (without the points at z = ±1) and Ti - T — T2- 

It is now easy to apply the method described in the previous section. At each time step the problem 

(2.4) is solved by replacing F 11 by F ( corresponding to F at t = n 5t ), T by T, and by taking 

. a 2 (2(3+1) a p 2 . . , 

A r = — + ^ r — — + ^ m (2.4e). 


dr 


dr 


It is easy to obtain T from (2.18), after having calculated T (except at r - 0). Note that the value of 
T at the axis can then be found by interpolation. 


2.3 Treatment of the singularity by omitting the axis 


In this section we propose another method to deal with the singularity which involves 
excluding the points at the axis. That can be done using the Gauss-Radau collocation points in the 
r-direction, between 0 and R, as suggested in [4, p.91], or by using an even number of Gauss- 
Lobatto collocation points in the r-direction, between -R and R, as proposed in [16] for a simple 
problem. Here we apply the latter technique to the general method described in section 2.1 by 
creating a new computational interval [-R,R], which corresponds to the diameter of the cylinder. 
Note that we are not forming the interval [— R,R] by extending each radial interval from [0,R] but 
by joining one radial interval [0,R] at the angle 0 with another at (0 + it ) for 0 6 [0 , tc]. Thus, 
the radial dependence of the functions is approximated by a Chebyshev expansion between — R and 
R on a well defined distribution of points which are physically meaningful. The singularity still 
remains when the axis belongs to the computational domain. However, by choosing an even 
number of Gauss-Lobatto-Chebyshev collocation points, the axis is automatically excluded. The 
problem of inforcing a boundary condition for the r-direction at the axis, when using Gauss- 
Lobatto points, is thus avoided. Furthermore, this device allows us to calculate the derivative with 
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*** 

respect to the r-direction in the right-hand-side F n of (2.3a) and the term *P in (2.4c) without 
taking into account the axis. For instance, for the evaluation of 'F we need to discretize the 
equation (2.4c). This requires the use of the formula (2.16) in which the collocation points, 
between -R and R, write as : 

ri = R^i , & = cos(jti / N) , i = 0,...,N. (2.22) 

It follows that in (2. 16a) : 


31 = 1 / R ( 2 - 23 ) 

From a practical point of view, values of N, such that N=5 n 3P (n and p being integers), must be 
chosen when using the FFT algorithm [14] ( to use the pseudospectral technique for the derivatives 
relative to the r-direction in F n ) in order to avoid the points at the axis. Of course, to remove this 
constraint on N, it is always possible to use formulas of type (2.11). In this case, the FFTs 
would be replaced by products of matrices (only for the calculation of the derivatives related to the 
r-direction). 

2.4 Assessment of these methods with an analytical solution 

In order to evaluate the above methods on problem (2.1)-(2.2), we have chosen the exact 


solution 


T ex (r,0,z) = F(r) H(0) G(z) + C , (Cel?), 

F(r) = r 7 (r -R) ; G(z) = (z 2 -1) (z 8 + 1) ; H(0) = sin 0 + cos(90 - 1) , 


(2.24a) 

(2.24b) 
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from which we define the Dirichlet boundary conditions 


Tip; = C , on Ti. 

(2.25) 

The velocities are then 


v Iex = sin0 r (R 2 - r) , 

(2.26a) 

v e ex = cos0 (2 r R 2 - r 2 (3 + 2z)) , 

(2.26b) 

v Zex = sin0 r (1 - z 2 ). 

(2.26c) 

It then follows that Ft is given by 


F t = V.VT - AT. 

(2.27) 

The chosen initial condition is 


T o (r,0,z) = T ex (r,0,z) + 80 (R - r)5 r (z 3 - z) sin 4 0. 

(2.28) 


The computations were carried out on a Cray XMP computer for the case C=1 and R=l. 

In the tables presented below, the residual of O (denoted R^}, the L 2 and Sup errors of 
(denoted E^ and E^ respectively) are displayed. Note that R^> j E^ 2 and E^ are calculated on the 

collocation points of and that ft = (T, v r , vq, v z ). Furthermore 



(2.29) 
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V N-l K M-l , , . xl9 

(N-l)K(M- l) ^ I l9 NKw( r i’ 0 k.Zj) - ^ex(ri»0k»Zj)| . ( 2 - 3 °) 

Efl P = Max i.kj | i^NKM{ri,9k,Zj) - 'Oex(ri> 0 k>zj)| , (2.31) 

with i=l,..,N-l , k=l,..,K , j=l,..,M-l ; $ ex = exact solution of and 

"Onkm = numerical solution of i!>. 

The problem (2.1)-(2.2) was solved by means of the following methods : 

- Method I described in section 2.1 by enforcing the real value of T at the axis (i.e. T=C at r=0), 

- Method II described in section 2.2, 

- Method m described in section 2.3. 

Note that method I is not applicable to real problems because the values of independent variables 
are not known at the axis. This method is nevertheless presented in this paper in order to assess the 
effect of the variable change used in the method II to treat the singularity (effect on the accuracy 
and rate of convergence). 

In all tables, we denote : 

- N r the number of collocation points between r = 0 and R, for methods I and II, 
and between r = -R and R, for method III. 

- K 0 the number of collocation points between 9=0 and 2k, for methods I and n, 
and between 0 = 0 and 0 = Jt for method III. 

- M^ the number of collocation points between z = -1 and 1 for all methods. 
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Table 1 presents results concerning the convergence (towards the steady state) and the 
accuracy of all three methods, after 200 time cycles. These results correspond to a parameter o 
(occurring in (2.3)) fixed to 0.5 and [3 equal to —1 (occurring in (2.18) and (2.19)). For 
approximately the same number of points, the accuracy is better with methods I and II than 
method in. With methods I and II, the machine accuracy (~ <9(10' 12 ) ) is reached with 
only 10x20x13 (i.e. N r x K e x M z ) points, after 200 iterations. For method in, when 
taking N r x K e x M z = 26 x 10 x 13, the error stagnates at C>(10' 8 ) after 200 time cycles. To 
decrease the error to 0(10'^) with this method, the number of points must be increased 
to N r x Kq x M z = 46 x 10 x 13. Despite the increase of points in the r-direction to N r = 46, the 
rate of convergence has not diminished significantly (residual R r r = 9.6 10"^ after 200 iterations). 
Whereas, in methods I and H, the increase of polynomials in the r-direction affects much more the 
speed of convergence. Table 1 shows that for these methods, the residual after 200 time cycles, 

is O(10' 4 ) when N r = 21. To reach an accuracy of <9(10‘ 12 ) with such number of points in the 

r-direction, 500 iterations are needed. When increasing the number of polynomials in the 
r-direction, the density of points near the origin, for methods I and n, is significantly increased as 
compared to the one in method IK. This leads to much stronger values of the term 1/r near the 
axis in the two first methods than in the last one. Therefore, for stability reasons, the time step for 
methods I and II is smaller than the one for method III. This restriction on the time step 5t 
decreases the speed of convergence of method I and II relative to method m. A rigorous stability 
study here is quite complex and would be beyond the scope of this paper. However, we must point 
out that the method III gives good results for solutions (2.24)-type for they do not involve strong 
gradients near the axis (in (2.24), dT/dr = 0 at the axis). When the solutions are stiffer close to the 
axis, method ffl is not as well adapted as method O, because the number of polynomials needed to 
represent the solution accurately becomes unacceptable. One possibility to overcome this difficulty, 
would be to use a mapping procedure as in [3,17], to concentrate the distribution of points in 
zones where sharp gradients in the solution are exhibited. For approximatively the same number of 
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polynomials, the CPU time per time step needed for technique III is about 1.5 times higher than 
the one required by the other methods. This difference is mostly due to the switches between the 
two kinds of coordinate system associated with this method (see section 2.3). These switches also 
generate programming difficulties and memory problems. Table 1 shows that methods I and II 
have the same order of accuracy but that the r ate of convergence is greater with II than I. This 
difference of speed is probably due to the better properties of the ADI operator (1 — o St Ar) 
following the variable change used in section 2.2. This difference has been found to be smaller 
when the number of points is higher. 

Finally, for more flexibility and efficiency, we have chosen the technique II to extend the 
method to Navier-Stokes and energy equations. By exploring different values of (3, we have noted 
that for P > — 1 the accuracy is poor and for P < — 1 the method is unstable (even when using very 
small time steps). In fact, the only value for which the method II works efficiently is P = -1. 
Table 2 shows the results of method n for different numbers of points and values of o. The value 
a = 0.5 was found to optimize the rate of convergence. In [3] for the solution of a diffusion 
equation in a domain between two concentric spheres the optimal value of o was 0.6. In the 
remainder of this paper we shall use P = -1 and o = 0.5. 

3. Numerical method for the Navier-Stokes and heat equations 
3.1 Description of the method 

In this section a method for solving Boussinesq-Navier-Stokes and energy type 
equations is presented. We use primitive variables [18] in a fixed reference frame. The problem to 
be solved will be the transport-diffusion type equation (2.1) coupled with 

= Pr Av z - V.Vv z - — + cos(atp + b) Pr Ra T + F Vl , 

9t dz 


(3.1a) 
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— = Pr Av r - V.Vv r - — + sin(ccp + d) sin(e0 + f) Pr Ra T + F Vr , 

dr 


(3.1b) 


^ V0 = Pr Av fl - V. Vv e - 1 — + sin(c(p + d) cos(e0) Pr Ra T + F V(J , 


at 


aa 


(3.1c) 


v. v =^ + ^ + i^ + ^=o , 

r dr r d0 dz 


(3. Id) 


where Av z and V.Vv z are obtained by replacing T by v z in (2.1b), (2.1c) and (2. Id), and 


2 -\ 2 _ 


dr 2 r dr r 2 r 2 dg 2 dz 2 r 2 d 0 


dvo 


(3.1e) 


\2 -\ 2 _ 


Avf 


= d v 9 t i dy e v e , 1 9 v e , d v e , 2 5v -~ 
3r 2 r dr r2 r 2 a0 2 az 2 r 2 ae 


(3.10 


__ 9v r Vfi 9v r v fl 2 9v r 

V.Vv r = v r — - + — — - - -®- + v z — , 

dr r d0 r dz 


( 3. 1 g) 


, r dvo Vo d Vo v r Vo 

V.Vvq = v r — 1 + -f— d - + -^r 9 - + v ; 


dr 1 d0 


dz 


(3 . 1 h) 


The above problem is solved in the domain Qi.with a, b, c, d, e, f, <p, Pr and Ra as parameters, 
and F Vl , F Vr , Fv e as given functions. The initial and boundary conditions will be defined later for 

each specific problem. The equations (2.1) and (3.1a,b,c) are solved by using the technique 
explained in section 2.2 with the variable changes (2.18) and (2.19) for T, v r , vq and v z , and 
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p(r,0,z) = rP p(r,0,z), when (3 = -1. For the problem (3.1), the typical numerical difficulties lie 
in the constraint V.V = 0 and the lack of boundary conditions for the pressure. In the steady case, 
these difficulties can be overcome by using the artificial compressibility method [2,5], which 
involves a perturbed continuity equation 

— + 6 V.V = 0 , (3.2) 

9t 

where 6 is a strictly positive constant. This equation has no physical meaning before the steady 
state 9 /3t = 0 is reached. Another possible technique to deal with the problem of pressure has been 
used in [19] to solve 3-D Navier-Stokes and energy equations in a parallelepiped domain. This 
method is analogous to the projection method [2] with an elegant iterative process. By applying the 
above variable changes, equation (2.1) transforms into (2.20) with T as unknown and the system 
(3.1a,b,c)-(3.2) becomes 

— = £ = r (Pr Av z - V.VvJ - ^ + cos(a<p + b) Pr Ra T + r F v , , (3.3a) 

at az 

— = G ■ r(Pr Av r - V.VvJ + £ * — , 

at ar 

+ sin(c(p + d) sin(e0 + f) Pr Ra T + r F Vf , (3.3b) 

^-=// = r(Pr Av 0 - V.Vv e ) - i-— + sin(c<p + d) cos(e0) Pr Rat T +rF v# . (3.3c) 

at r 30 


1^ + £V.V = 0 , 

r Ot 


(3.3d) 
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with the diffusive term rAv z and the convective term rV.Vv z obtained by replacing T by v z in 
(2.20b) and (2.20c), and 


r Av r = 




r 0r r 2 de 2 


| c?Vt 2 3^9 
+ 9z 2 r2 30 ’ 


(3.3e) 


r Av 0 = 


0 VQ i 9vq 

~dS~ T dr 


O V0 


ar 


l a v 9 | 2 dv r 
3z 2 r2 ao ’ 


(3.3f) 


r V.Vv r = y 


0 Vr 

0r 


Vr 2 ! V9 V9 2 , V z 3 y r > 
r 2 r 2 dd r2 r dz 


(3.3g) 


r V.Vv e = 


v r dve 
r 0r 


, ve 3ve 
r 2 de 


+ 


v z dv e 

r a z ’ 


(3.3h) 


V.V 


l^L + X§Xi + i^ 
r ar r 2 90 r 9z 


(3.3i) 


The new system to be solved is equation (2.20) coupled with equations (3.3a,b,c,d), with 
T, v r , ve, v z and p as unknowns. The system will be solved by the following steps : 

(i) First, we calculate p n+1 by using the following explicit scheme on (3.3d) : 


p" +1 = p" - 5t £ r (V.V) n 


(3.4) 
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(ii) With T n , v?, Vq, and v5, we deduce F" from (2.20) and then find T n+1 , by replacing 

F n by r, T by T , and by taking A r = — - ^ in (2.4). 

dr ™ r 

_ i , ~n,n+l ~n,n+l ~n,n+l 

(iii) Using T n+1 , p n+1 , v?, v£, and v^, we then calculate E ,G and H 
relative to (3.3a), (3.3b) and (3.3c) respectively. 

(iv) We then deduce v* +1 by replacing F n by F n ’" +1 , T by v z , and use the same operator 
A r as for T n+1 in (2.4). 

(v) For the calculation of v? +1 we replace F n by G , T by v r and set 

in(2 - 4) - 

dr 

(vi) Finally, Vq +1 is computed by replacing F n by H , T by vq and by taking the 
same operator A r as for v? +1 in (2.4). 

The right-hand-side F n corresponds to F at the time t = n St. The computation of 
E n ' n+ \ G n ’" +1 and // n ’ n+1 is made by using the velocities at the time t = n 5t, and the pressure 
and the temperature at the time t = (n+l)5t. For the stages (iv), (v) and (vi), the operators Ae, A z 
and A r are multiplied by Pr. 

3.2 Evaluation of the method with an analytical solution 

In order to assess the method described in Section 3.1 on problem (2.1)-(3.1), we 
have chosen the exact solutions (2.24) and (2.26a,b,c) for T, v r , vq and v z , respectively, and 


p ex = cos(20) r 3 (1 - z 2 ). 


(3.5) 
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From these exact solutions, we readily deduce the Dirichlet boundary conditions of the variables 
and the forcing terms Ft (see 2.27), F v ,, F Ve and F Vr . The chosen initial conditions are (2.28), for 

T, (3.5) for p, (2.26a, c) perturbed with 

P v = 20 (R - r) (z 2 - 1) cos 59 , ( 3 - 6 ) 

for v r and v z , and (2.26b) perturbed with r Pv for v e . The calculations have been done with 
Pr = Ra = C = R=l,a = b = c= e = 0, f = d = 7t/2. 

Table 3 shows the convergence and accuracy results of the method, after 2000 time 
cycles, when 8t = 10' 2 and N r x Kg x M z = 10 x 20 x 1 1. After this number of iterations the 
maximum of accuracy has been almost reached. In these calculations we noticed that the maximum 
error was always located at the collocation points close to the axis. This is certainly a consequence 
of the variable change used to remove the singularity. We remark that the rate of convergence is 
connected to the arbitrary constant £ (occurring in (3.3d)). For this specific problem, with this 
number of polynomials, and this value of 8t, the value E = 30 was found optimal. An interesting 
possibility for improving the convergence, would be to employ a technique analogous to [4, 
p.149] to adjust the coefficient £ at each time step. Note that the initial conditions are totally 
different from the final steady solutions. In real physical problems we generally choose the initial 
conditions to be as close as possible to the expected ones. 

3.3 Evaluation of the method on two simple physical problems 

In this section we shall apply the above method to two simple convective flow problems 
of incompressible fluids using the real Navier-Stokes and heat equations within the Boussinesq 
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approximation. The goal is not to study these physical problems but to assess the method by 
comparing it with two other ones. Thus, we shall have to solve the equation (2.1) with Ft = 0 
(corresponding to the heat equation) coupled with equations (3.1) when 
p v = p v ^ = F Vi = b = d = f = 0 and a = c = e = 1 (corresponding to the Navier-Stokes equation), in 

the domain Sip These equations are written in dimensionless variables with the Prandtl number 
Pr and the Rayleigh number Ra as physical dimensionless parameters. For the comparison 
with other methods, we have chosen two physical problems which have 2-D solutions. The 
first generates an axisymmetric solution and the second one a solution independent of the 
axial direction. The other methods chosen for the comparison, were designed for solving 2-D 
Navier-Stokes and heat equations. The Rayleigh number Ra has deliberately been chosen to be 
small, in order to avoid problems related to the accuracy of each method. In forthcoming sections 
we shall denote the maximum values of v r , vq, v z , v x and v y by v rmax , vq max , v zmax , 
v x max 111(1 v y max- respectively. The variables v x and v y are the velocities in Cartesian 
coordinates such that v x = v r cos 6 — vq sin 0 and v y = v r sin 0 + vq cos 0. As these maxima 
are calculated at the collocation points particular to each method, small differences are expected in 
the results. 

3.3.1 Axisymmetric solution 

For this problem, the fluid is differentially heated in a vertical cylinder having rigid walls 
as boundaries. The axisymmetric convective flow is generated by an axisymmetric non-linear 
temperature distribution at the initial state. The equations described earlier were solved with <p =0 
(i.e. gravity parallel to the axis r = 0), Pr = 1, and the boundary and initial conditions 


v r = ve = v z = 0 on T] , 

T = 0 at z = 1 , T=1 at z = -l , T = (z- l) 2 /4 at r = R , 


(3.7a) 

(3.7b) 
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v r = VQ = v z = 0 at t = 0 , (3.8a) 

T = (z- 1) 2 /4 at t = 0. (3.8b) 

Table 4 gives the values of v r max and v z max , at different Rayleigh numbers, obtained with the 
finite volume program PHOENICS [20] in cylindrical coordinates, the finite element code FTDAP 
[9,10] in Cartesian coordinates and our present method. The calculations were done by taking 
N r xM z = 18x36 with PHOENICS, N r xM z = 17x31 with FIDAP and N r xM z = 13 x 19 
with our method. Note that for PHOENICS and FIDAP the mesh was regularly spaced. We can 
observe a perfect agreement between our results and those given by FIDAP and a small 
discrepancy with PHOENICS which yields about 5% higher values. The error associated with the 
first-order finite volume approximation used in [20] has the effect of slightly increasing the values 
of v r max and v z max . Note that these results were obtained with an optimum mesh size when 
using FIDAP and PHOENICS. Fig.2 exhibits the velocity and temperature fields when Ra = 150. 
We can see that v z max is located on the axis. 

3.3.2 Solution independent of the axial direction 

In this section, the fluid is heated from the side in a horizontal cylinder with rigid walls. 
The convective flow is buoyancy generated by taking cp = 7 t / 2 (i.e. gravity perpendicular to the 
axis r = 0), Pr = 1, and the boundary and initial conditions 


v r = V0 = v z = 0 

at the sidewalls (r=R) , 

(3.9a) 

3v r _ 3 v 9 _ 3v z 

= 0 at the endwalls (z = + 1) , 

(3.9b) 

dz dz dz 



T = cos 0 

on Ti , 

(3.9c) 
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v r = vg = v z = 0 at t = 0 , (3. 10a) 

T = rcos0 at t = 0. (3.10b) 

By enforcing v z = 0 in the equations, the solution becomes independent of the axial direction. 
Table 5 shows the values of v x max and v y max , obtained with Ra = 60 and 120, first by using 
FDDAP and then our method. The calculations with FIDAP were done by taking a regularly spaced 
Cartesian grid 29 x 29 and N r x Kg = 13 x 18 with our method. As in the previous section, we 
obtain a good agreement between these two methods. Fig.3 presents the velocity and temperature 
fields when Ra = 120. 

4. Application to a natural convection problem 

In this section the method is illustrated with a convection problem of a differentially 
heated fluid in a horizontal cylinder having a length 2H and a diameter D (A=D/2H is the aspect 
ratio). The vertical endwalls of the cylinder are maintained at constant temperatures Ti and T 2 
while the horizontal boundary is conducting. The convective flow is generated by the buoyancy 
force as soon as T^Ti. The intensity of the resulting flow is connected to the magnitude of the 
difference T 2 -T 1 . All boundaries of the cylinder are rigid with no-slip conditions for the velocities. 
The mathematical model is given by the same Navier-Stokes and energy equations as used in 
section 3.3 where <p = rc / 2 (i.e. gravity perpendicular to the axis r = 0). These equations are cast 
in dimensionless form with H, H^/k, k/H, and T 2 —T 1 as characteristic scales for lenght, time, 
velocity and temperature, respectively. The resulting dimensionless computational domain is 
with R = A. The physical dimensionless parameters are the Prandtl number Pr = v / k and 
the Rayleigh number Ra = gaH 3 (T 2 -Ti) / kv where g is the gravitational acceleration, a the 
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thermal expansion coefficient, v the kinematic viscosity, k the thermal diffusivity. The 
dimensionless boundary and initial conditions are 


v r = vq = v z = 0 at the endwalls (z = ± 1) and the sidewalls (r - R) , (4.1a) 

T = 0 at z = 1 , T=1 at z = -1 , T = 0.5 (1 + z) at r = R , (4.1b) 

v r = vq = v z = 0 at t = 0 , (4.2a) 

T = 0.5 (1 + z) at t = 0. (4.2b) 


Table 6 shows the values of v xmax , v ymax and v zmax , obtained with Pr= 1, A = 1 and 
Ra = 60, 100 and 600, first by using the finite element code FIDAP and then our method. For 
FIDAP we have employed an irregularly spaced Cartesian grid with N x x N y x N z = 17 x 17 x 1 1, 
similar to the example depicted by a figure in vol.3, section 35-27, of the FIDAP manual [9]. We 
denote N x , N y and N z the number of nodes in the x, y and z directions respectively. For our 
method all calculations (tables 6 and 7) have been performed with a spatial resolution defined 
by N r x K 0 x M z = 13 x 18 x 19 and a time-step 5t varying between 3. 10' 3 and 1. 10‘ 3 
according to the stiffness of the solution. The product 8t £ occurring in (3.4) is maintained at the 
value 0.2 for all computations related to this section. Table 6 shows good agreement between these 
two methods. Table 7 summarizes the details of our computations corresponding to Pr = 1, A = 1 
and Ra between 1000 and 13000. The starting condition (4.2) has been used only for the lowest 
Rayleigh number. For the others, the starting condition is the solution calculated with a lower 
Ra. With the chosen number of polynomials, on a Cray XMP computer, the computational cost for 
one time-step is 0.081 second. Since the goal was to illustrate the method and its applications 
rather than undertake an exhaustive study, we did not try to explore the solutions for 
Ra > 13000. The physical study of this problem has been carried out in [21-23] with other 
parameters by means of a finite difference method. Table 7 shows in particular that the time step 
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decreases with increasing Ra which leads to a corresponding increase of the number of cycles. 
However, even with 6000 cycles the total CPU time remains reasonable (486 sec.). Note that for 
high Ra we have stopped our calculations when the divergence (corresponding to (3. Id)) 
V.V * 10' 5 . For these values of V.V, we noticed that the values of the velocity and temperature 
fields do not change significantly when decreasing the divergence to lower values. This fact was 
independently pointed out in [24]. Figures 4 and 5 exhibit the velocity and temperature fields when 
Ra = 13000. The projection of these fields in the (r.z)-planes is shown in figure 4 for 0 equal to 
0°, 40° and 80°. At 9 = 0°, the flow has one roll in each comer. At 40° and 80°, it consists of 
roughly one cell with a contraction in the middle. Fig. 5 exhibits the fields at various cross- 
sections (r,0). At z = 0.76, i.e. close to the hot endwall, the flow is ascending and perfectly 
symmetrical to the descending flow at z = —0.76, i.e. close to the cold endwall. At z = 0, it is 
mainly directed from the side to the middle with two secondary flows in the vicinity of the wall 
crossing the median axis parallel to the gravity. Note that for Ra = 13000, the flow encounted 
could be defined as a " Boundary Layer Driven Regime ", as proposed by [23]. 

5. Conclusion 

A Fourier-Chebyshev pseudospectral method has been proposed for solving steady 3D 
Navier-Stokes and heat equations in cylindrical cavities. The principal computational difficulties in 
such geometries come from the presence of the singularity at the axis and the number of collocation 
points req uir ed to accurately represent the solution. Concerning the first difficulty, two techniques 
have been considered and evaluated. One is based on a variable change and the other on a special 
distribution of the collocation points. The technique using the variable change has been found more 
adapted to the general method. Moreover, the use of a pseudospectral method admits an accurate 
representation of the solution with a small number of collocation points. 
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This paper has been focused on the solution of steady problems. Therefore, difficulties 
connected to the incompressibility condition and the pressure could be efficiently circumvented by 
means of the artificial compressibility method. The latter, coupled with the use of the ADI 
procedure in the general method leads to a very fast and accurate algorithm for computing steady 
3-D incompressible flows. However, for an extension of this method to unsteady problems, the 
present way to deal with pressure could be time consuming, due to the resulting iterative process. 
A better technique to extend it to time dependent problems, might be the use of a semi-direct 
method as in [19], to satisfy the divergence free field. Finally, the application of our method to a 
natural convection problem, has shown that this technique can cope efficiently with relatively large 
Rayleigh numbers. 
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Table headings 


Table 1 - Convergence and accuracy after 200 time cycles for methods I, II and HI. 


Table 2 - Convergence and accuracy for method II. 


Table 3 - Convergence and accuracy for problem of section 3.2. 


Table 4 - Accuracy for problem of section 3.3.1. 


Table 5 - Accuracy for problem of section 3.3.2. 


Table 6 - Accuracy for problem of section 4. 


Table 7 - Results regarding the problem of section 4. 


Table 1 




10x20x 13 9.0 10-3 7.65 lO' 10 1.24 10' 11 8.36 10' 11 

I 21x20x 13 4.0 10-3 1.22 10' 4 1.33 10' 6 2.47 10'5 


10x20x13 1.2 1 O' 2 4.96 10" 12 8.38 10* 13 7.37 10’ 12 

II 21 x20x 13 4.0 10-3 1.05 10- 5 7.32 10' 7 2.00 1 O' 6 

26x 10x 13 1.6 10- 2 1.82 10" 11 2.92 1 O' 9 1.72 10‘ 8 

HI 46x 10x 13 1.0 10- 2 9.60 lO' 10 2.51 10' 11 1.80 lO’ 10 
















\ 


Table 2 


N r x Kg x M z 

a 

8t 

Number 

of 

cycles 

Rj 

Ej 


9 x 16 x 9 

0.5 

2.0 lO- 2 

50 

8.21 10-7 

3.98 10-3 

1.76 lO' 2 

10 x 20 x 11 

0.5 

1.2 lO' 2 

200 

4.68 10- 12 

4.12 lO’ 13 

1.23 10' 12 

10 x 20 x 11 

0.75 

1.2 lO' 2 

300 

1.60 10' 11 

1.16 10' 11 

1.07 lO’ 10 

10x20x 11 

1.0 

1.2 lO' 2 

400 

7.24 10’ 9 

9.37 lO' 9 

1.01 lO' 7 

13 x 24 x 13 

0.5 

0.8 lO' 2 

300 

1.78 10- 12 

1.68 10- 12 

5.64 10' 12 











Table 3 


e Rf) E^ z 

30 9.13 10-13 7.63 10' 13 2.23 10‘ 12 

25 9.92 10' 13 7.74 10' 13 2.29 10' 12 

20 1.12 10' 12 8.33 lO ’ 13 2.43 10' 12 

30 2.19 10-9 1.47 10-8 i .53 10‘ 7 

25 7.52 lO ’ 9 5.85 10’ 8 6.30 10' 7 

20 2.38 lO ' 8 2.24 10’ 7 2.48 10' 6 

30 5.47 10- 11 2.46 lO ' 10 1.57 10' 9 

25 1.87 10-10 9.51 lO'lO 6.29 10’ 9 

20 5.91 10 - 1 ° 3.52 10' 9 2.42 10' 8 

30 9.08 10 -H 8.70 10 - 1 ° 6.42 10' 9 

25 3.06 10-10 3.37 10' 9 2.59 10' 8 

20 9.56 10 - 1 ° 1.25 10' 8 1 . 00 10' 7 

30 6.04 10-6 1.32 lO * 7 9.86 10' 7 

25 2.07 10* 5 5.19 10‘ 7 4.00 10‘ 6 

20 6.59 lO ’ 3 1.94 10-6 1.56 10' 5 




















Table 4 


Ra 


PHOENICS 


Vz max 


v rmax 


FIDAP 


v z max 


v r max 


SPECTRAL 


Vz max 


Vrmax 


10 

100 

150 


1.77 10- 2 
1.91 10’ 1 
3.02 10’ 1 


4.15 10' 3 
4.50 10-2 
6.90 IQ ' 2 


1.66 10-2 
1.80 lO ’ 1 
2.82 10' 1 


3.98 lO ' 3 
4.28 lO ' 2 
6.64 lO * 2 


1.66 lO ' 2 
1.80 10* 1 
2.82 10- 1 


3.96 lO ' 3 
4.23 lO ' 2 
6.59 lO ’ 2 
















Table 5 


fidap spectral 

Ra v x max Vy max v x max v y max 

60 1.39 1.41 1.38 1.40 


120 2.66 2.69 2.62 2.68 









FIDAP 


SPECTRAL 


v x max 


v y max 


v z max 


Vx max 


v y max 


v zmax 











Table 7 


10000 
11 
12 
13 


8t 

Number 

of 

cycles 

V.V 
(3. Id) 

v r max 

v 0 max 

v z max 

3.0 10- 3 


1.1 10- 4 

7.82 

7.91 

9.89 

3.0 10-3 


5.0 lO- 5 

11.72 

12.45 

14.82 

3.0 10-3 

4000 

2.2 10-5 

14.29 

16.20 

18.26 

3.0 10-3 

4000 

3.0 10-5 

17.05 

19.45 

20.99 

3.0 10-3 

4000 

4.0 10-5 

19.40 

22.09 

23.29 

3.0 10-3 

4000 

5.0 10-5 

21.53 

24.80 

25.35 

2.0 10-3 

4000 

6.4 10-5 

23.44 

27.33 

27.39 

2.0 10-3 

4000 

8.8 10-5 

25.18 

29.59 

29.26 

2.0 10-3 

4000 

1.3 lO- 4 

26.77 

31.62 

31.01 

1.5 10-3 

5000 

3.7 lO' 5 

28.24 

33.44 

32.65 

1.5 10-3 

5000 

6.1 lO’ 5 

29.62 

35.61 

34.21 

1.5 10-3 

5000 

9.6 10*5 

30.90 

37.69 

35.68 

1.0 10-3 

6000 

2.9 10-5 

32.15 

39.66 

37.09 


Initial 

condition 


Sol. (4.2) 

Sol. at Gr=1000 
Sol. at Gr=2000 
Sol. at Gr=3000 
Sol. at Gr=4000 
Sol. at Gr=5000 
Sol. at Gr=6000 
Sol. at Gr=7000 
Sol. at Gr=8000 
Sol. at Gr=9000 
Sol. at Gr= 10000 
Sol. atGr=11000 
Sol. at Gr= 12000 












Figure legends 


Fig. 1 Geometry of the enclosure and frame of reference. 

Fig. 2 Problem of section 3.3. 1 : (a) velocity field, (b) temperature field for Ra = 1 50. 

Fig. 3 Problem of section 3.3.2 : (a) velocity field, (b) temperature field for Ra = 120. 

Fig. 4 Problem of section 4 : (a) velocity field, (b) temperature field, in the (r,z) planes, 

at Ra = 13000 for 9 = 0°, 40° and 80°. 

Fig. 5 Problem of section 4 : (a) velocity field, (b) temperature field, in the (r,0) planes, 
at Ra = 13000 for z = 0.76, 0 and -0.76. 
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